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ABSTRACT: This work presents a methodological approach where distribution of potential pollution sources and 
historical frequency of pollution incidents are included in the assessment of risk to groundwater resources. Calibrated flow 
and transport models were setup, and integrated with a FORTRAN coded risk model. The risk model was run to generate 
synthetic contaminant source terms that were subsequently transported by the coupled flow and transport models. The 
ranges of calibrated values obtained for horizontal and vertical hydraulic conductivities are 5.787x1 0' 6 - 2.31 5x1 0' 5 m/s and 
5.787x1 0' 8 - 1.1 57x1 0' 7 m/s, respectively. Layering within the aquifer is thought to account for the large vertical anisotropy. 
The corresponding values for the specific yield and specific storage are 0.10 - 0.12, and lxlO' 4 - 5xl0' 4 , respectively. The 
estimated historical frequency of pollution occurrence per stress period of 91 days is 35, while the probability of pollution 
occurring at any of the sources over a period of stress period is calculated as 0.0014. The generated risk maps provide 
anticipatory capability for risk quantification and show that approximately 55% of the monitoring boreholes have high (i.e. 
70 - 100 %) likelihood of having contaminant concentration being less than 0.06 mg/l throughout the simulation. 
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I. INTRODUCTION 

The statistics of groundwater use (Zektste and Everett, 2004) has shown that groundwater has become one of the most 
important global natural resources. In Denmark, Malta and Saudi Arabia, groundwater constitutes the only major source of 
water supply, while in Tunisia and Belgium, approximately 95 and 83 % respectively of the country's water resources is 
sourced from groundwater. This importance is underpinned by the inherent valuable properties of the resource compared to 
the surface water. These properties include higher quality, better protection from anthropogenic activities, less seasonal 
variation, large storage, relatively inexpensive developmental cost and wide spread occurrence, among others. However, 
occurrences of global pollution is increasingly becoming one of the major environmental concerns that pose threats to 
groundwater resources and can potentially cause depletion of these values. Cases of pollution occurrences have been 
reported in Beijing in China (Marilyn, 2001), Taejon area in Korea (Jeong, 2001), as well as in the developing countries 
(Egbu, 2004). The scale of impacted groundwater bodies is becoming more widespread and the persistence of groundwater 
pollutants is increasing, making the cost of aquifer restoration to be excessive in many cases. Therefore, pollution 
prevention rather than pollution management appears to proffer the basis for a truly sustainable approach to groundwater 
protection. 

Generally, the extent of protection of groundwater resources is assessed based on the risk posed by anthropogenic 
activities, or the measure of ease with which an infiltrating pollutant reaches the groundwater resource. The contemporary 
techniques for the assessment of risk and vulnerability to groundwater can be broadly be categorized into three groups 
namely: ranking index methods, process-based computer simulation, and post-pollution assessment methods (Al-Adamat et 
al, 2003; Aller et al, 1985; Connell and Dale, 2003; Foster, 1987; Gustafson, 1989; Van Stempvoort et.al, 1992; Khan and 
Liang, 1989; Rao et al, 1985; and USEPA, 1989). According to Haimes (2006), risk is defined as the result of a threat with 
adverse effects to a vulnerable system. Following from this definition, it implies that the most vulnerable groundwater is not 
at risk without the presence of threats. This concept of risk therefore comprises of two-dimensional components. From the 
point of view of the water resources, the first component is the probability of contaminant source terms being generated from 
potential pollution sources at the ground surface while the second component is the assessment of impact of such occurrence. 
The ranking index methods are essentially qualitative and subjective in approach, and generally lacking good scientific 
judgement. Also, none of the existing methods incorporates quantifications of the probability of occurrence of polluting 
source terms from the potential sources. 

In this work, potential pollution source is defined as a single, identifiable and localized source with risk of 
discharging pollutants that may infiltrate into the aquifer. Hence, the fact that potential sources often constitute threats to 
groundwater resources necessitates the need to include their occurrence and distribution in quantifying any risk posed to the 
underlying groundwater resources. Quantifying the probability of occurrence of pollution is generally the most uncertain part 
of the assessment of risk posed to groundwater system. Therefore, this work is set to demonstrate the field application of a 
novel two-dimensional risk assessment methodological approach, where assessment of risk to groundwater resource 
incorporates both the quantification of the probability of occurrence of source terms, as well as the impact of such pollution 
event. A more detailed discussion on the concept and structure of the method has already been presented in Oladeji and Elgy 
(2012). 
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II. METHODOLOGY 

The algorithm of the methodology implemented in this work is presented in Figure 1 . The risk model components 
of the Risk Assessment Method (RAM) structure are implemented by a computer program written in a standard FORTRAN 
90. The risk model is run over the same period of time as the flow and transport models. Contaminant source terms are 
generated for each stress period of the simulation, which are then transported within the subsurface environment using the 
coupled flow and transport models. The effects of the generated source terms are assessed by observing the spatial and 
temporal concentrations of the contaminant at pre-determined monitoring boreholes within the aquifer, as well as by 
counting the number of times the contaminant concentration exceeds user defined ranges of concentration magnitude. Two 
grid systems namely a local system and a global system are implemented in the risk model. This is to allow the risk model to 
be run independently of the flow and transport models, as well as to increase the efficiency in the implementation of the risk 
model in terms of the time and memory requirements. The local grid may be equal to or smaller than the global grid system 
used in groundwater flow and contaminant transport models. 

The purpose of the Risk Model Parametization module (see Figure 1) is to prepare the input data for the risk model. 
This involves obtaining locations of the potential pollution sources at the ground surface. The geographical locations of the 
potential sources are initially obtained using global positioning equipment, and subsequently transformed into the 
appropriate layer, row and column numbers of the model grid using GIS utilities. 
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Figure 1: Conceptual algorithm of the risk assessment method 

The process of the generation of synthetic source terms by the risk model involves determination of the condition 
for the generation of a synthetic source term, as well as the magnitude of the associated contaminant mass. In order to 
determine the condition for the generation of a synthetic source term in any stress period, the risk model computes a 
parameter called Probability of Pollution Occurrence (PPO) for each of the active model cells containing potential pollution 
source(s), using Equation 1 : 

No of days of previously occurred pollution events/ 

„ „ „ / No of sources 

PPO = - (1) 

Total no. of days under consideration 



Note that equation 1 is implemented when the PPO has to be evaluated during the execution of the risk model. If PPO is 
already known by priori knowledge, then the values can be entered appropriately. 

Next, for each day in the current stress period and for each active node, the risk model generates a Random Number 
(RN) between 0.0 and 1.0, and compares the random generated number (RN) with the PPO value computed. If PPO > RN, 
then a pollution incident is assumed to occur. However, an array of integer values of between 0 and 5 is set for each model 
cell, and these dynamically control the actual proportion of infiltrating contaminant mass per each event. For example, a 
value of 0 indicates that all the synthetically generated probable mass loading rate of contaminant at a particular active node 
are transported from that source, while a value of 5 indicates that, although synthetic contaminant mass is generated at the 
source, no pollutant is assumed to be transported from that source. This utility is used to incorporate source control 
capability as mitigating measures based on the ground conditions. Also, in order to determine a contaminant mass associated 
with each of the synthetic pollution event, the risk model samples a random value from a range of the minimum and 
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maximum values of contaminant mass loading rates that have previously occurred based on historical records or other priori 
knowledge. 

In this work, the probability distribution of the historical pollution events is represented by a Poisson distribution. 
A Poisson distribution expresses the probability of a number of events occurring in a fixed period of time if these events 
occur with a known average rate and independent of the time since the last event. If the expected number of occurrence in 
an interval of time is given as A, then the probability that there are exactly k occurrences (where k being a non-negative 
integer) is given as: 



where e is the base of the natural logarithm, and k! is the factorial of k. 

Finally, the risk model extracts the spatial and temporal contaminant source terms generated for all the stress 
periods, and presents the source terms in a format that is compatible with the transport model input data file. In addition, the 
risk model transforms the row and column numbers from the local grid system into that of the global grid system used for 
flow/transport model. 

The generated source terms serve as input into the transport model. After each run, the risk model counts the 
number of times at which the contaminant concentrations exceeds ranges of user defined concentration magnitudes, and 
express the result as a risk. These steps represents the first realisation of the risk assessment method, and subsequently 
iterated such that the frequency distribution of the observed contaminant concentration at points of interest in the aquifer is 
developed by observing the number of times that the user defined concentration magnitude is exceeded. 

2.1 Description of the study area 

The study area (Figure 2) extends 17 km from south to north, and 13 km from west to east, covering approximately 
221 km 2 , and forms part of the larger West Midlands conurbation. Birmingham has been a major manufacturing centre for 
over a century, and contains many industrial processes that may constitute risk to groundwater. According to Powell et al 
(2000), the area is underlain by the Triassic sandstone, and consists of fluvioglacial thick successions of the sandstone 
deposit, commonly referred to as the Sherwood sandstone Group. The basal aquifer unit is the Kidderminster sandstone 
Formation, overlain by the Wildmoor and Broms grove sandstone Formations. These three sandstone aquifer horizons have 
distinct hydrogeological properties. The superficial geology consists of till, sand, gravel and alluvium, and overlies the 
Triassic geological successions. The soil types present within the study area are highly varied and complex, and their 
distribution and composition are influenced by factors such as climate, geology, geomorphology and hydrology. The 
landuse pattern is largely urban. The major rivers present within the study area are shown in Figure 2. 

A major geological structure within the study area is the Birmingham Fault, which juxtaposes the Triassic mudstone 
to the east against the Triassic sandstone to the west. According to Allen et. al, (1997), transmissivity is generally reduced 
across the Birmingham Fault, and Knipe et. al, (1993) modelled the zone across the fault as a reduced aquifer thickness in 
order to represent the reduced transmissivity across the fault. The Triassic sandstone Group form the major aquifers in the 
study area, and their hydrogeological characteristics are dominantly controlled by lithological variations, vertical 
heterogeneity, anisotropy, fractures, scale of measurement, as well as uneven variations in the aquifer thickness (Allen et al, 
1997). The inorganic, organic contamination and the occurrence of pollution related acidification within the urban aquifer of 
Birmingham area have been respectively studied by Ford and Tellam, (1994), Rivett et al (1990), and Ford et al (1992). 

2.2 Numerical groundwater flow and transport model 

The U.S. Geological Survey numerical finite-difference groundwater flow model MODFLOW 2000 (Harbaugh et 
al. 2000) forms the basis for the model development and the subsequent calibration process. The initial conditions for the 
model setup are presented in Table 1. For the purpose of model calibration, the model was set up to run under transient 
conditions covering 20 years from January 1970 to December 1989, and validated using groundwater head data spanning 
from March 1985 to February 2015. The eastern and the southern boundaries of the aquifer geometry were defined using no 
flow cells because of the presence of Birmingham Fault which is assumed to inhibit flow across it. Also, the western 
boundary was defined as no flow because of the presence of the Westphalian Formations, which are older non permeable 
geological material. The northern boundary was delineated by the groundwater divide along an anticlinal axis on the base of 
Triassic sandstones (Knipe et al, 1993), and therefore also represented as no flow boundary. 

The initial groundwater heads within the aquifer horizons were interpolated from the available groundwater level 
observation data. Where the river path exists, the initial groundwater head was defined by the river stage. The rivers present 
within the model domain (see Figure 2) are contained within the first model layer. The Modflow head dependent river 
package allows the model to calculate the amount of streambed percolation and groundwater inflow to each river reach, 
using river bed conductance and the difference between the calculated model cell hydraulic head and the stage of the river. 
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Figure 2: Location of the study area 
Table 1 : Summary of initial and final flow model data 



Input parameter 


Initial Conditions 


Boundary Conditions 


No flow conditions for the west, 
east, north and south boundaries. 


Spatial discretization (m) 


No of rows: 760; No of columns: 
600; Ax=25; Ay=25 


Horizontal hydraulic 
conductivity (m/s) 


Broms grove Fm =5.0xl0~ 5 
Wildmoor sst Fm =2.0x10" 
Kidderminster Fm =3.0xl0~ 4 


Vertical hydraulic 
conductivity (m/s) 


Broms grove Fm =5.0x10"' 
Wildmoor sst Fm =1.0xl0" 6 
Kidderminster Fm =5.0xl0" 7 


Specific yield 


Bromsgrove Fm = 0.10 
Wildmoor sst Fm = 0.10 
Kidderminster Fm = 0.10 


Specific storage 


Bromsgrove Fm =lxlO" ;i 
Wildmoor sst Fm =5xl0" 3 
Kidderminster Fm= lxlO" 3 


Effective Porosity (%) 


Layer 1 = 26.8 %; 
Layer 2 = 23.8 %; 
Layer 3 = 26.4 % 


rj L Or ; a V T / Ovl; Mol. 
Diff. Coeff. 


20 m; 0.1 m; 0.01; 0.0 m 2 /s 



The river bed conductance controls the rate of flow to or from the river according to the difference between the 
river stage and the modelled groundwater head in the uppermost active cell. River bed conductance also defines the 
maximum leakage rate to the aquifer when heads fall below the bed of the river. The river paths contained within the model 
boundary were divided into reaches. The interaction with the underlying aquifer was simulated between each reach and the 
model cell that contains that reach. The river bed elevations were estimated from the digital elevation map (DEM), and the 
river stage was assumed to be 0.5 m above the estimated river bed. 

The constant flux well package was used to simulate pumping from the 12 abstraction boreholes (gw_l - gw_12) 
located within the model domain (Figure 2). The average abstraction rates used for each stress period are presented in 
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Figure 3a. The dearth of data for the abstraction rates restricts the use of the actual values of the historical abstraction rates in 
the transient model calibration. Also, the three locations (sw_l - sw_3) where surface water abstraction is being abstracted 
within the study area (Figure 2) are lined and therefore no interaction with the underlying aquifer is assumed. The daily 
recharge values were estimated for the study area using a recharge calculator based on series of spreadsheet calculations 
underpinned by the Food and Agricultural Organisation (FAO) methodology. The average annual recharge values for the 
study area is 112 mm/yr, and represents the final recharge value distributed across each stress period (Figure 3b) in the flow 
model using the Modflow recharge package. 

The model was calibrated with hydraulic head data spanning over 20 year period, from March 1970 and February 
1989, by minimizing the residuals between the observed and the simulated groundwater head data. The convergence 
criterion for the hydraulic heads was set to 0.01 m within the Preconditioned Conjugate Gradient 2 (PCG2) solver package 
The PCG2 solver package iteratively refines the initial estimates of the model until either an acceptable measure of residual 
is achieved or the difference between the results of successive iterations is less than a user specified convergence criteria. 
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(c) GW head against the residual after the model calibration (d) Source terms generated during the simulation 

Figure 3: Abstraction and recharge as well as the residual and source terms 

The groundwater level monitoring data obtained from the observation boreholes were used as targets for assessing 
the effectiveness of the calibration process. In order to make future predictions of the risk posed to the groundwater 
resource, the flow model was validated and further set up to run from March 1985 to February 2015. The corresponding 
groundwater head computed at the beginning of March 1985, which corresponds to the elapsed time of 15 years of the 
calibrated flow model, was set to be the initial groundwater heads for the model layers under the predictive transient 
simulation. Also, a predictive transport model (Zheng and Wang, 1999) was setup under transient conditions covering 30 
years from March 1985 to February 2015. 

2.3 Risk model 

The records of the historical occurrence of chloride related pollution incidents within the study area spanning over 
January 2002 - December 2009 were used to determine the pattern of the historical occurrence of pollution incidents. This 
data is considered to be sufficient in the demonstration of the applicability of this risk assessment method, but thought that 
longer record length will be more representative for the study area. In order to keep the model simple, this work considered 
pollution sources from only five common industries within the Birmingham area, and these include food, chemical, garages, 
mining and mineral industries. The analysis of the historical pollution data shows that the average number of pollution 
incidents from the selected industries for the study area was 35 events per stress period of 91 days. These values were used 
to compute the PPO for the study area. The total number of pollution sources within the study area was 28 1 . Using equation 
1, the probability of pollution occurring at any of the sources over a given stress period of 91 days is computed as: 
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No of days of previously occurred pollutionevents / 35/ 

PPO = /Noof sources = Jm = QQQU 

Total no. of days under consideration 9 1 



The calculated PPO was applied to each of the active cell of the risk model for each stress period. Although, the method as 
demonstrated in this work assumed that each type of the industry has the same probability of occurrence of pollution 
incidents in each of the stress period, this method can be adapted to assess the risk from the individual type of industry with 
varying probability of pollution occurrences per stress period for each of the industry. 



III. PRESENTATION AND DISCUSSION OF THE RESULTS 

The results of the applications of the risk assessment method are discussed under the following headings. 

3.1 Flow model 

The model was setup as a three-layer model, which represents (from top) the Bromsgrove, Wildmoor and the 
Kidderminster sandstone Formations. The respective final horizontal hydraulic conductivity values (K h ) are 5.787xl0" 6 , 
2.315xl0~ 5 , and 3.472xl0" 5 m/s. The corresponding values for vertical hydraulic conductivity (K v ) are 5.787xl0~ 8 , 1.157x10" 
7 , and 5.787xl0~ 8 m/s, respectively. Hydraulic layering within the Permo-Triassic sandstone aquifer is thought to account for 
the large vertical anisotropy in the final calibrated hydraulic conductivity values. The final horizontal conductivity values 
are within the same range compared to the values obtained by Allen et al. (1997) from field test data, for the respective 
aquifer horizons, as well as to those values obtained by Knipe et al. (1993). Furthermore, the final values for the specific 
yield are 0.12, 0.10 and 0.12, and for the specific storage are lxlO" 4 , 5xl0~ 4 , and lxlO" 4 , respectively for Bromsgrove, 
Wildmoor and Kiderminster Formations. The value reported by Allen et al. (1997) for the specific yield of the undivided 
Sherwood sandstone Group is 0.12. Knipe et al. (1993) and Rushton and Salmon (1993) respectively reported specific yield 
value of 0.15, and 0.10 for the Bromsgrove sandstone Formation, and these values are similar to those obtained in this work. 

The plot of the observations and the corresponding residual for the model calibration is presented in Figure 3c. The 
percentage numerical error associated with the volumetric balance is less than 0.1 % throughout the duration of the 
simulation, and this indicates that the model was converging without any significant numerical error. A sufficient degree of 
match was obtained between the measured head observations and the simulated equivalents (Figure 3c). The flow model 
was subsequently setup as a predictive model for a 30-year period spanning March 1985 - February 2015. The March 1985 
groundwater head output of the calibrated model was used as the initial water level for the predictive model. 

3.2 Risk model 

The risk model was set up over the same period of 30 years, from March 1985 to February 2015. The global grid 
system for the risk model is the same as that of the flow and transport model. The number of potential sources of chloride 
pollution within the case study area was scoped for the food processing, garages, quarries, mineral and chemical industries, 
and 281 locations were identified as sources of chloride pollutant. The average occurrence of chloride pollution incidents 
over the eight-year period (2002 - 2009) was estimated as 35 incidents per stress period (91 days). Hence, the probability of 
pollution occurrence for each of the potential pollution source was computed to be 0.0014. The range of the contaminant 
mass loading rate of the historical pollution incidents within the study area is estimated based on the available information to 
be 250 - 5000 kg per incident. This represents the range of values within which the contaminant mass loading rate was 
sampled on each occasion when a synthetic pollution incident occurred. 

The risk model was run to generate source terms for the transport model. The number of pollution events generated 
per each day is summarised in Figure 3d, and it indicates that zero incident occurred in 9,453 days out of the total 10,920 
days. Single daily events occurred in 1363 days, two daily events in 98 days, while three daily events occurred in only 6 
days. Given the generating mechanism used to create pollution events, the frequency of occurrence of 0, 1, 2... events per 
day should follow a Poisson distribution. The Chi (X2) test verifies this, as it showed that the probability that these are by 
chance from different distributions is very remote (1.02937xl0~ 6 ), and therefore confirms that the mechanism for generating 
random events is correct. 



3.3 Transport model 

The predictive transport model was setup over a 30-year period (March 1985 to February 2015). Chloride 
contaminant was used to demonstrate the applicability of this method because of its conservative and non reactive nature 
within the natural subsurface environment. The initial background chloride concentration was set to zero in order to prevent 
the background concentration from masking the effects of the generated source terms. The mass balance error is generally 
less than 0.1 %. The spatial distribution of the contaminant after 30 years of simulation is presented in Figure 4, and shows 
that the contaminant concentration within the study area varies between 0.0 - 0.1 mg/1. High correlation exists between the 
simulated contaminant concentration (see Figure 4) and the pattern of the distribution of the potential pollution sources (see 
Figure 2). 

3.4 Risk assessment results 

The risk to groundwater resources within the study area was determined as the number of times which chloride 
concentrations at the 12 observation boreholes exceeded the user defined concentration intervals. The two user defined 
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concentration intervals used are < 0.06 and > 0.06 mg/1. The risk of the contaminant concentration not exceeding 0.06 mg/1 
within the first model layer is presented in Figure 4. That is, Figure 4 shows that there is 70 - 100 % likelihood that the 
contaminant concentrations within 55 % of the monitoring boreholes will be less than 0.06 mg/1 throughout the simulation 
period, while the confidence of this likelihood is less than 70 % at all other monitoring points. A major limitation in the 
application of this method is that it has not been possible to calibrate and validate the probabilistic component of the risk 
model. Also, it is very difficult to determine the probability of occurrence of extreme events because of lack of records. 
Furthermore, high pollution incidents usually have higher correlation compared to low pollution events, and the capability 
for assessing this type of correlation is not incorporated in the current level of the development of this risk assessment 
method. 



402000 41 0000 




402000 41 0000 



Figure 4: Final contaminant concentration and risk map for 0.06 mg/1 

However, this approach of risk assessment presents a more resolved definition of risk distribution within the study 
area. The existing groundwater vulnerability maps classified the study area as high vulnerability area. Considering the 
potential pollution sources present within the study area, this work has further resolved this generic classification, with 
capability to delineate both spatially and temporally, the areas that are prone to quantified risk from the probable pollution. 
This involves isolating the specific risk to the individual water source as well as the relative trend of the risk distribution 
within the study area. The roles played by the incorporation of the occurrence and distribution of potential pollution sources 
at ground surface are not accounted for in the contemporary risk assessment methods. 

Although, it is generally difficult to make definitive statements about the predictive accuracy of one risk assessment 
method compared to other, this method does offer a greater insight in the quantification of risks to a groundwater source, and 
will be preferred in circumstances where the anticipatory assessment of the risk to groundwater resources is required. 
Generally, any risk assessment approach requires continuous and dynamic reviewing and update, and this process may span 
over several years. Also, it is generally difficult to test regional risk assessment methods on a field scale in a similar way as 
site specific risk assessment methods will be tested or validated. Alternative validation approach is to apply more than one 
method to the same site and compare the results. 

IV. CONCLUSIONS 

This work presents the demonstration of the application of a two dimensional risk assessment method using a case 
study. The method assesses the effects of a probabilistic source terms, based on the historical frequency or externally 
determined probability distributions of pollution incidents, on groundwater source. Integrated flow, risk and transport 
models were setup and used as predictive tool to assess effects of synthetically generated source terms at the observation 
boreholes. Risk is calculated as the number of times at which a user defined contaminant concentration magnitudes are 
exceeded over a period of time. The risk of pollution from a number of sources all occurring by chance together was 
evaluated, and this capability to combine the risk to a groundwater feature from numerous potential sources of pollution 
proved to be a great asset to the method. 
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